Holstein polaron in the presence of disorder 
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Non-local, inhomogeneous and retarded response observed in experiments is reproduced by in- 
troducing the Inhomogeneous Momentum Average (IMA) method to study single polaron problems 
with disorder in the on-site potential and/or spatial variations of the electron-phonon couplings 
and/or phonon frequencies. We show that the electron-phonon coupling gives rise to an additional 
inhomogeneous, strongly retarded potential, which makes instant approximations questionable. The 
accuracy of IMA is demonstrated by comparison with results from the approximation free Diagram- 
matic Monte Carlo (DMC) method. Its simplicity allows for easy study of many problems that 
were previously unaccessible. As an example, we show how inhomogeneities in the electron-phonon 
coupling lead to nonlocal, retarded response in scanning tunneling microscopy (STM) images. 



PACS numbers: 71.38.-k, 72.10.Di, 63.20.kd 

Understanding the nature of the materials under the 
focus of current basic research, as well as the develop- 
ment of novel applications, is unambiguously linked to 
the physics of quasiparticles in disordered systems and 
coupled to bosonic modes of different origins. Thus, the 
manganites which exhibit colossal magnetoresistance are 
doped materials [T with considerable electron-phonon 
(el-ph) as well as electron-magnon and electron-orbitron 
couplings [2 ; the interplay between disorder [3] and the 
coupling to bosons manifests itself in the peculiarities 
of their phase diagram [4 . Similarly, the underdoped 
high temperature cuprate superconductors are inhomo- 
geneous materials [5^ with rather strong [6 and inhomo- 
geneous [7 coupling to phonons. Another example con- 
cerns charge transport in organic thin- film transistors [8] , 
which is dominated by polaron jumps between different 
potential traps [9 . The importance of the interplay be- 
tween disorder and coupling to boson fields is magnified 
by the fact that weak coupling to a boson mode that is 
relatively unimportant in a clean system may result in 
dramatic effects in disordered compounds [10 . 

The tremendous difficulties in treating the problem of a 
polaron in the presence of even a single impurity resulted 
in the 30 year delay between the first results based on the 
adiabatic approximation [11 and the recent approxima- 
tion free solution by the DMC method [10]. However, 
the current level of technology requires theoretical pre- 
dictions for numerous systems whose inhomogeneity is 
not limited to a single impurity but includes any form 
of spatial inhomogeneity, plus potential wells or barriers 
representing surfaces, interfaces or quantum dot wells, 
besides cases where both the energy of the bosonic mode 
and its coupling to the electron are inhomogeneous [12 . 
Solving such general problems for large systems by nu- 
merical methods is still effectively impossible. 

In this Letter we study accurately yet efficiently the 



single Holstein polaron problem with disorder in the po- 
tential, the strength of the coupling constant and the 
frequency of the phonon by developing the Inhomoge- 
neous Momentum Average (IMA) method. The method 
is based on the Momentum Average (MA) approximation 
used to study translationally invariant systems, like Hol- 
stein [13 and more general models [TT. IMA takes any 
potential inhomogeneity into account exactly and can also 
handle spatial variations of the coupling constant and of 
the frequency of the boson modes. Comparing results 
of IMA with approximation free data from DMC allows 
us to gauge its accuracy. We perform this comparison 
in one-dimension (ID) where the worst accuracy of IMA 
is expected [13 . The IMA approximation can then be 
systematically improved [15 so that in combination with 
DMC fW" one gets an accurate and fast tool to study all 
the systems described above within the framework of a 
controllable and efficient approximation scheme. 

Given the low computational cost of IMA, a rapid scan 
of large regions of the parameter space is now possible. 
To compare with experimental observations, we compute 
STM images of inhomogeneous systems of large enough 
sizes to render DMC studies impractical, due to the enor- 
mous computational cost required for the analytic con- 
tinuation method [16 . This allows us to prove the nonlo- 
cal nature of a system's response to inhomogeneities and 
demonstrate the importance of the retardation effects. 

To avoid cumbersome expressions we consider one im- 
purity in an otherwise homogeneous ID system, indicat- 
ing generalizations where suitable. The Hamiltonian is: 

n = no^Uo^ T/ei-ph (1) 

where Ho = + ^'^') + ^Zli ^l^i ^^^^ 

part, Uq = —UcqCq the on-site attraction to the impurity 
placed at site z = 0, and T4i-ph = ^Xli^l^i (^^1 + 
the el-ph interaction. The electron's spin is irrelevant. 
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The goal is to compute the retarded Green's function: 
nr ■ ^ /nl n, ^ tin\ V- (0|c.|a)(a|cj|0) 

a 

where G{u;) = [co — H -\- ir]]~'^ and h = 1^ because it has 
information on the single electron eigenstates 7i\a) = 
Ea\<^)- Also, the local density of states (LDOS) A{i^u;) = 
— ;^ImG(z, i, cj) is measured directly by STM. 

We first introduce two additional Green's functions. 
One is the free electron Green's function 



Go{iJ,u;) = {0\c,Go{u;)c]\0) 



1 

2^ 



dk- 



^ik{Ri — Rj) 



e/c + 



where = —2tcosk for nearest neighbor hopping. Gen- 
eralizations to other dispersions and higher dimension- 
ality are trivial. The second is the "disorder" Green's 
function Gd(^, j, Cc;) = (0|QGd(^)c]|0), corresponding to 

= ^\g=o = Ho-\-Uo' Using Dyson's identity Gd{uj) = 
Go{uj) + Gd{u;)UoGo{u;) straightforwardly leads to: 



Gd{iJ,uj) = Go{i-j,uj) - U 



Go{i,Lu)Go{j,u) 
l^UGo{0,u) ' 



(3) 



since Go{i^j^u;) = Go{i — j^uo) = Go(j — h^) (the sec- 
ond equality holds if time reversal symmetry is obeyed). 
Eq. ([3| is valid in any dimension, if the appropriate Gq 
is used. For more complicated disorder potentials one 
can find Gd(^, by a suitable generalization: Hd is a 
quadratic Hamiltonian and can always be diagonalized. 
Thus, Gd is known and treats the on-site disorder exactly. 

Since Ti = Hd + Ki-ph, we can now proceed to calcu- 
late the desired Green's function in terms of Gd(^, j, c^;). 
Using Dyson's identity once, we find: 

G{iJ,uj) = Gd{iJ,uj)^g^Fi{iJi,ij)Gd{jiJ,uj) (4) 



With this approximation only the functions survive 
in Eq. ([5|, whose general solution is then Fn{i,j^u;) = 
An{j,uj)Fn-i{i^ j^uj) [13]. The continued fractions: 



ngGd{jJ,uj- nft) 



1 - gGd (j, j, ^ - nn)An^i (j, uj) 



(6) 



are simple to compute. We now insert cj) = 

74i(j, cj)G(z, j, cj) in Eq. Q, resulting in G(z,j, cj) = 
Gd(^,J» + ^Eji ji,cj)Ai(ji,cj)Gd(ji, To 
make this system of coupled equations converge fast with 
the cutoff in ji, we define 



An{u;) = An{j,u;)\u^o = An{j,u;)\ 



(7) 



since if = 0, then Gd{jJ,u;) = Go(j, = Go(0,cj) 
irrespective of j. The same holds for finite U but \j\ 
oo, since sites located very far from the impurity are not 
sensitive to its presence at the origin. Introducing the 
"effective interaction" potential: 



^0 ( j, ^) = gAi ( j, uj) - Sma(o) (^) 



(8) 



where = {0\c,G{u)c]b]''\0), with Fo(i,j» = 

G{i^j^uo). Using the Dyson identity again, we find that: 

F„{i,3,uj) =gJ2 Gd(ii,i,w - nn){0\ciG{io)c]^b]^b]''\0) 

^gGdUJ^^- nQ) [nFn-i{iJ,uj) + . (5) 

Within the IMA^^) approximation, we set in all these 
equations Gd(ji, j, — n^l) ^ if j 7^ ji and n > 0. This 
is a good low-energy approximation, because the ground- 
state (GS) of the polaron is below the spectrum of Hd and 
so for uj ~ Eqs^ Gd{i^j^^ — nQ) decreases exponentially 
with the distance \i — j|, the decrease being steeper for 
larger n. In other words, IMA^^^ ignores exponentially 
small terms. Like MA^^\ IMA^^^ is in fact accurate at 
all energies because it obeys multiple sum rules [131 US] • 
It also becomes exact both for ^ ^ and t ^ 0. 



where the bulk MA^^^ self-energy is T^-^j^io) {co) = gAi{uj) 
[15], we can rewrite the equation for G{i^j^u;) as: 

G{iJ,u;) = Gd(i, j,a;)+^G(i, ji,cj)i;o(ji,cj)Gd(ji, 
31 

(9) 

where uj = u — Sma(o) (^)- This equation is very efficient 
to solve numerically, because vo{j^u;) rapidly with 
increasing In fact, a cutoff \j\ < 5 suffices for con- 
vergence, although a cutoff of or 1 does not, showing 
that vo{j^uj) is spread over a few sites around the impu- 
rity. Note that inhomogeneities in the values of g and 
Q are easy to deal with, since one simply has to use the 
appropriate gj and Qj values in Eqs. (|6| and (|8|. 

Equation ^ reveals a two-fold role of the el-ph inter- 
action. If the solution was just G{i^j^uo) = Gd(^, 
it would mean that the renormalized quasiparticle - the 
polaron - interacts with the bare impurity potential Uq. 
However, the second term shows that the impurity po- 
tential is renormalized as well 



Uo Uo ^^vo{j,uj)c]cj, 



(10) 



and is no longer local and has significant retardation ef- 
fects through its cj-dependence. In other words, because 
of el-ph interactions, the dressed polaron interacts with 
a renormalized, retarded disorder potential. 

Since Uq is treated exactly, we expect the validity of 
this approximation to mirror that of the MA^^^ for pure 
Holstein (17 = 0) therefore to worsen in lower dimensions, 
when g ^ t (such that the effective ID el-ph coupling 
A = g'^ / {2tQ) ~ 1; remember that the approximation be- 
comes exact for both A ^ 0, 00) and for smaller ft [T5] . 
Even then, a larger U improves the accuracy as it pushes 
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FIG. 1: (color online) (a) and (c) Ground-state energies; 
and (b),(d)-(f) Spectral weights at the impurity site vs. the 
impurity potential [/, for t = 1, Q = 2 and g = 1.5 so 
that A = 0.5626, in (a) and (b), respectively Q = 0.2 and 
g = 7^2, ^/O^S, VOA so that A = 0.5, 0.75, 1 in (c)-(f). 



the OS to lower energies. However, for small U and Q and 
for A ~ 1 we need to go to a higher level of the approxi- 
mation. Like for the homogeneous MA(i) solution [15], m 
we neglect the exponentially small contributions 
only if there are n > 2 phonons present. The IMA^^^ 
solution is similar to Eq. (|9|, except Co is now renor- 
malized by the bulk Sma(i) (^) self-energy [15 , while the 
renormalized potential is also more accurate: vi{j^uj) = 
g'^Xj^^/ [1 - gxj^^[A2{j,uj) - Ai{j,ij - n)]] - Sma(i)(^) 
where Xj^^j = ^ma(o) (j? J? oo — ^) is the MA^^^ solution 
of Eq. ([9| at a shifted energy. If necessary, one can go 
to a higher IMA^"^^ (n > 2) level in the same way. 

We gauge the accuracy of IMA for model ([T]) against 
DMC results. Fig. [T] shows the OS energy and quasipar- 
ticle weight Zgs at the impurity site vs. the impurity po- 
tential U. The agreement is very good in (a) and (b) even 
though these are ID results, where IMA is least accurate. 
This is partially due to the large ft = 2t used [15 . For a 
worst-case scenario, we plot in (c)-(f) results for a much 
smaller ft/t = 0.2, for weak, medium and strong cou- 
plings. Now we see clear differences, although they are 
quantitative, not qualitative. For £^gs, the disagreements 
at small U just mirrors the errors of MA for the Holstein 
model [15]. As expected, as U increases and Eg^ moves to 
lower energies, the accuracy improves. Zgg shows more 
significant errors at small [/, with IMA overestimating 
the correct answer. This is not surprising, since for such 
small Q one expects many phonons to be created at many 
sites and a higher level n of IMA is needed. However, 
even levels n = 0, 1 capture the physics quite well. More- 
over, given the spectral weight sum rules satisfied exactly 
(6 for IMA(O), 8 for IMA^^), see Ref. [15 ), we expect the 
spectral weight at all energies to be similarly accurate. 

In Fig. [l] (a), (b) we also show what happens if we 
set the additional potential vo{j^u;) 0, i.e. we use an 




FIG. 2: (color online) (a) ID LDOS A(i,E). The polaron 
band and part of the second bound state band are shown. 
The vertical dashed lines mark the edges of the disordered 
region. The slanted lines show the dispersion of some features 
in the spectrum; (b) Analog of (a) for a 2D sample. The 
plot shows A(i,j = 25,0;) vs i; (c,d,e) show the 2D A(i,j,E) 
vs. i,j for E = 0,0.25 and 0.5 respectively. In all plots the 
energy is measured from the ground-state polaron energy of 
the corresponding uniform system with A = 1, Q = 0.5t. 



"instantaneous" approximation (see e.g. Refs. [rzj.fTS]) 
where the el-ph coupling is assumed to renormalize the 

2 J. 

potential Uq ^ Uq — ^"^jCjCj. This overall shift is 
given, in IMA, by the bulk self-energy Sj^^co) ^ ^gs) ^ 
—g'^/Q through the renormalized Co. Clearly, this works 
quite poorly even for this rather small A, and becomes 
considerably worse as g and therefore vo{j^uj) increase. 

In Fig. [2] we show IMA results for a nontrivial prob- 
lem whose treatment by the DMC method is not fea- 
sible due to enormous computational costs. Here we 
study the site and energy dependent map of the LDOS 
of an area with randomly chosen el-ph effective coupling 
Xi = gf /{2tQ) G [0.9, 1.1] embedded in an otherwise uni- 
form infinite system with A = 1. We do not add an on-site 
disorder potential, although it can be included exactly as 
discussed above. The LDOS maps in Fig.[2];a),(b) show 
clear evidence of the nonlocal response of the system to 
such inhomogeneities in the coupling constant, with var- 
ious features changing their position in space as the en- 
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FIG. 3: (color online) (a) Fourier transform of the 2D LDOS 
map at energy E = 0.06t above the bulk ground-state; (b) 
Same for ky = and E/t = 0.03,0.06,0.12. The lines show 
predicted scattering kg vectors. For more details, see text. 



ergy is changed and '^0(^7 ^) varies (the slanted lines show 
some examples). This is more pronounced in ID, since 
in 2D some of these shifts proceed in the direction per- 
pendicular to the line with j = 25. The retarded nature 
is apparent not only through the different-looking LDOS 
maps at different energies inside the disordered region, 
but also by the slow convergence towards the uniform 
bulk value and the Friedel-like oscillations seen in the uni- 
form region surrounding it. As expected, the wavelength 
of these oscillations decreases with increasing energy. 

The 2D LDOS maps at different energies show no cor- 
relations, but their Fourier transforms (FT) should show 
peaks at k values where scattering between qp states of 
that energy is likeliest. Such analysis is well known for 
STM data in cuprates, for scattering both on impuri- 
ties |19| and inhomogeneities in the el-ph coupling [7 . 
In Fig. [3]; a) we show the FT at £; = 0.06t above the 
bulk polaron ground-state (we averaged over FT of 200 
LDOS maps such as shown in Fig. [2j but for 50x50 site 
inhomogeneous regions. We removed the k = peak for 
clarity). At such low energies, the bulk polaron disper- 
sion E{k) = -2r(cosfc^ + cosky - 2) ^ n^\s?/(2m*), so 

-k 
L86 

for A = 1, 17 = 0.5t. We show ks as a dashed line found to 
be in agreement with the FT data, as confirmed in Fig. 
[sj^b) for several energies. This demonstrates that it is the 
polaron and not the bare particle that is scattered by in- 
homogeneities. We believe that the second peak visible 
for higher energies is due to inelastic scattering between 
the first and second polaron bound states, but this needs 
further study. In any event, these results show agreement 
with the general phenomenology seen in the experimen- 
tal data [7, 19^ well beyond the single site disorder case 
usually considered theoretically [20^. 

Because IMA is very fast (2D LDOS maps such as 
shown in Fig. [2|c) take about 50s to generate on an or- 
dinary desktop; moreover this type of calculation is ideal 



we expect signal up to a /cg = 2y 2m*E /fi^ for k 
scattering. The effective polaron mass is m* /m 



to parallelize, with different CPUs for different energies) 
one can easily study very large disordered regions, for dif- 
ferent types of disorder in the on-site potential (whether 
short- or long-range) and/or el-ph coupling or phonon 
frequency, in any dimension. Such studies will reveal the 
effects of each kind of disorder and the way in which mul- 
tiple types of inhomogeneities do or do not "interfere". 
Moreover, such studies can be extended to other bare qp 
dispersions as well as models beyond Holstein, such as 
systems with multiple phonon modes or with electron- 
phonon coupling dependent on the phonon momentum, 
where the bulk MA solutions are already known [14 . 

In conclusion, we have studied the Holstein polaron 
problem with disorder by developing an accurate, con- 
trollable and fast computational method suitable for a 
huge class of problems that were unaccessible until now. 
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